Possibilities of the Discrete Fourier Transform for Determining the 
Order-Chaos Transition in a Dynamical System 

Carlos R. Fadragas(l), Ruben Orozco-Morales(2) , Juan V. Lorenzo- Ginori(2), 
(l)Physics Department, Central University of Las Villas, Santa Clara, Cuba, 
(2) Centre for Studying Electronics and Information Technologies, 
Central University of Las Villas, Santa Clara, Cuba. 

Abstract 

This paper is devoted to a discussion of the Discrete Fourier Transform (DFT) 
representation of a chaotic finite-duration sequence. This representation has 
the advantage that is itself a finite-duration sequence corresponding to samples 
equally spaced in the frequency domain. The Fast Fourier Transform (FFT) 
algoritm allows us an effective computation, and it can be applied to a relatively 
short time series. DFT representation requirements were analized and applied 
for determining the order-chaos transition in a nonlinear system described by 
the equation x[n + 1] = ra;[n](l — x[n]). Its effectiveness was demonstrated by 
comparing the results with those obtained by calculating the largest Lyapounov 
exponent for the time series set, obtained from the logistic equation. 

1 Introduction 

One of the most studied problems in the scientific research deals with the pro- 
cessing of a time series a;[2], . . . , x[N], composed by an experimental data 
sequence, spaced in time domain, sometimes evenly, sometimes unevenly. Such 
a sequence is obtained by sucessively sampling over a dynamical observable, 
characterizing a dynamical system under investigation. The goal is, as a result 
of processing, to obtain information about the behavior of the dynamical system 
that produced the time series. When we deal with a linear dynamical system, 
obtainnig the information is relatively easy. However, there exist some difficul- 
ties when the dynamical system is a nonlinear one. In this case, the behavior 
of the dynamical system is rich enough. This was known for a long time, but it 
was not until the latest 20 years that the development of modern technologies 
and the digital numerical computation allowed finding solutions to important 
problems of the nonlinear dynamics, and to develop new methods in this field. 
The nonlinear dynamical analysis of a time series is a modern tool for scientific 
research, that has demonstrated its usefulness when it is to be applied in many 
branches of human knowledge. 

The analysis procedure begins by reconstructing a phase space portrait of 
the attractor by using delayed coordinates where the dynamics of the system 
can be observed|^, ^ |^, |l^, |2^. For doing this, vectorial representations using 
delayed coordinates from a scalar time series x[n]; n — 1,2, . . . , N , (N being the 
length of experimental data sequence), are made up. The embedding dimension 
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in the phase space, De, determines the dimensionahty of the vectorial field X[k\ 



X[k] = {x[k],x[k + T\,x[k + 2t], . . . ,x[k + {De - l)r]) (1) 

It is usual that the time delay, r, is an integer number times the signal sampling 
period. The maximum number of state vectors that can be reconstructed with 
delayed coordinates, given a experimental data sequence of length N, equals 
M ~ N — {De — 1)T. The vector set X[k]; k — 1, 2, . . . , Af , each with coordinates 
x[k + iT];i — 0, 1, 2, . . . , De — 1, signals to points on the orbit of the system 
attractor in the phase space, so that the time evolution of the system in the 
phase space configuration can be observed by taking the initial condition X[l] as 
a starting point and going on the behaviour of the vector state X[k]. This phase- 
space portrait is a representation of the behaviour of the dynamical system and 
can give us some information about the dynamics of the system[^, [n], 1^, |, 



0, 1^. The selection of the values of reconstruction parameters, and r, can 
be made by using various methods (| 0, ||, |, |, § |, 

After determining the reconstruction parameter values, one proceeds to cal- 
culate dimensions, exponents, and entropies. However, these quantities are 
difficult to compute, require a large data set, and degrade rapidly with additive 
noise. Generally, algorithms assume the repeatition of a same task for defining 
an adecuate value of the parameter that is being investigated Many 
authors recognize that the application of metric tools of the non linear dynam- 
ical analysis implies a high computational cost. It must be noticed that there 
are certain requirements related to the minimal length N of the time series to 
be processed in order for the results obtained by applying nonlinear dynamical 
analysis to be realiable|8|. 

Besides, according to the bibliography the transition from a regular to a 
chaotic behaviour has been studied in many dynamical systems, including the 
biological type ones. One of these is a fiuid fiow, where by calculating nonlinear 
dynamical parameters is possible to determine the order-chaos transition. CM 
Van den Bleek y JC SchoutenQ determined the order-chaos transition in a flu- 
idized bed by calculating correlation dimension, D2 and Kolmogorov entropy K 
as a function of the Reynolds number and the superficial gas velocity, indicating 
two different regimes, order and chaos, and the intermediate region. 

On the other hand, the discrte-time Fourier transform (DFT) representation 
and the digital technique for its processing, known as fast Fourier transform 
(FFT) has been widely used for analysing a time sequence and a linear time- 
invariant system. However, it is required to distinguish three types of time series: 
deterministic, chaotic, and random. This classification is based upon the value 
that the Kolmogorov entropy takes in each case|4j. The DFT representation of 
a deterministic time series has been widely used, and the DFT representation of 
a random time series there doesnt exist. We are interested now in analysing the 
possibility of applying the DFT representation to a chaotic finite-duration time 
series. We assume that the behaviour of a nonlinear dynamical system which 
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produces a time sequence n — 1, 2, 3, . . . , iV, can change from a order to 
a chaotic regime and viceversa, and that is reflected in a any way in the time 
series observed Q. 

It is postulated in this work that every method, that allows us to give an 
evaluation about the behaviour of a dynamical system, and that can act over a 
relatively short time series, has a great practical importance. 

The aim of this research is to evaluate the ability of the DFT representation 
of a discrete-time series for obtainning some criteria about the ordered-chaotic 
transition related to the behaviour of a dynamical system, without applying a 
priori discussed and laborious nonlinear dynamical analysis algorithms, char- 
acterized all of them by a high computational cost. In particular, the method 
proposed in this paper will be realiable if it is applied to discrete-time series set 
generated by a given dynamical system, where the comparative criteria are to 
be the goal. It is important to say that this method does not pretend to susti- 
tute those of the nonlinear dynamical analysis, but it can be used for analysis 
complementation. So that, the DFT representation can be applied to a shorter 
time series than those required by the nonlinear dynamical analysing. 

We produced a 120 time series set with 1024 data points in each time series, 
and for each of them we obtained its 1024 FFT representation. By observing 
the amplitude spectrum and the phase spectrum behaviour with the control 
parameter of the system, a well defined value of the control parameter (r = 3.57) 
can be determined, beyong which the behaviour either the amplitude spectrum 
and the phase spectrum became quite irregularly, in contrast with the regular 
behaviour of these magnitudes in the region before the indicated value above. 
There exists a complete correpondence when this result is compared with that 
obtained from the bifurcation diagram and from the largest Lyapounov exponent 
graphic, both as a function of the control parameter r. 

2 Methods and Materials 

2.1 Obtention of the discrete-time series set 

The logistic equation is a one-dimensional mapping of the real axis in the interval 
(0,1), that is a prototype of a nonlinear dynamical system widely used|^ and 
can be formulate as 

x[n + l]=rx[n]{l-x[n]), (2) 

and with the initial condition x[l] especified, one can obtained a time sequence 
of a given length by a recursive action. The character of the behaviour of a 
data sequence obtained with the iterative process can be modified by selecting 
conveniently a value of the control parameter r. In this work a value of the initial 
condition a;[l] — 0.65 was taken, and the interval r £ [2.8,4.0] was selected. A 
set of evenly spaced values of r, with a step S = 0.01, was taken to produce 
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a family of 120 discrete-time series, each containning N — 2^° samples. The 
building up of the histogram for each time series allows us to give inmediately 
a simple statistical characterization of these time series. 

In order to do much easier further computations data points were organized 
into a rectangular matrix containing 120 columns each of which is a time series 
with the length N = 1024 data points. 



2.2 Determining of the largest Lyapounov exponent Ai 

Lyapounov exponents spectrum deals with average exponential divergence or 
convergence of two neighboor orbits in the phase space. It is common to or- 
der the spectrum from largest to shortest exponents: Ai, A2, A3, . . .. A system 
containning one or more possitive exponents is to be defined as chaotic. For 
calculating the Lyapounov exponents stectrum some algorithms have been de- 
veloped, according to the particular situation, and being one of the most used 
algorithm that reported by Wolf et al[|l]. From the exponents spectrum, the 
largest exponent, Ai, decides the behaviour of the dynamical system. For calcu- 
lating the largest exponent one can refer too the algorithm proposed by Rosen- 
stein et al [pT| . For estimating Ai one can monitor the large term evolution of 
two neighboor orbits on the attractor reconstructed. Details of these algorithm 
can be found in refered papers. 

We determined the largest Lyapounov exponent, Ai, for each time series in 
the set, using a professional software ijisl. The selection of the time delay, r, and 
the embedding dimension, D,,, required for reconstructing the phase space is 
part of the problem and it is necessary to consider some selecting criteria. This 
selection can be made by different methods [|[ |l^, [l^, ^,J|, ^ |l^. For select- 
ing r one can cite: (a) Optimal filling of the phase space[g, (b) by determining 
the position of the first local minimum of the data autocorrelation function||l|, 
(c) by taking the first minimum in the mean mutual information plot pT[ , (d) 
based on the optimal tradeoff between redundance and the irrelevance pW, and 
by analysing of the differential equations of the system if they are given ]2^. 
For selecting can be indicated: (a) by applying the Grassberger-Procaccia 
method, which allows one to obtain simultaneously both correlation dimen- 
sion and embedding dimension p^, (b) by applying the false nearest neighboors 
method [T^, [sj, and (c) by applying geometric considerations [p^ . 

A time series from the logistic equation was analyzed by using the method (a) 
for T determining and by taking a value De = 3, which satisfies the Grassberger- 
Procaccia algoritm condition De > 2_D2, and Takens theory |l^, being D2 



0.5 a value reported by Hoyer et al|14 for this time series type. Figure 4 shows 
the effect of optimal filling of the phase space for one of the time series obtained 
here. 



Literature refers] 13 19|] that selecting embedding dimension as 

De>D2 + 1, (3) 
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is suffient to have a good representation for the attractor in the phase space, 
reconstructed in delayed coordinates from a scalar time series observed. Any 
reconstruction made with a value of the embedding dimension below the minimal 
value will produce a projection, in that dimension, of the original attractor, and 
therefore it will be more difficult for the interpretation of that projection or it 
will produce false results when a estimation of nonlinear dynamical parameter 
is to be made. It is important to take a value for embedding dimension close to 
2D2 + 1 because of it is not the case, computational calculation will be excesively 
made. 

2.3 Application of DFT representation: Amplitude spec- 
trum and phase spectrum. 

The DFT representation of a finite-length discrete-time series x[ri\;n = 0, 1, 2, ... , N— 
1, is given by 

N-l 

= nT. m^'^^^"'^''^; n = 0, 1, 2, . . . , iV - 1 (4) 
fe=o 

N-l 

X[k\= ^ x[n]e-^'(2'^/^)'=";A; = 0,l,2,...,iV-l (5) 

n=-0 

The sequence X[k];k = 0,1,2, ... ,N — 1 can be efficiently computed by a dig- 
ital algorithm known as Fast Fourier Transform (FFT), from which the am- 
plitude spectrum, \X[k]\ or abs{X[k]), and the phase spectrum, arg{X[k]) or 
angle{X[k]), are both obtained. 

Next wc obtain the DFT representation for each time scries in the set by 
the 1024 points-FFT of the set matrix. Since each time series produced by the 
logistic equation (with r e [2.8,4.0]) is a finite-duration bounded sequence, it 
satisfies the existence requirement for the DFT representation, that is 

00 

E kWI<oo (6) 

n=— 00 

We considered that it is convenient, for some purposes, to eliminate the dc 
component of signals produced by the logistic equation, before the amplitude 
spectrum mesh plot is made. 

2.4 Results 

Figures, from 1 to 3, show a sample of three representative time series and their 

corresponding histograms. In order to give a better view of each time series it 
was only considered 128 data points of each discrete-time series in its plotting 
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but it was completely considered when its corresponding histogram was made. 
In either case, the particular value of r parameter was indicated. 

The two-dimensional plot of the rectangular matrix produces a figure known 
as Feigenbaum or bifurcation diagram. It reflects the transition of a dynamical 
system towards chaos by a period-doubling process. Feigenbaum diagram are 
obtained when parameter values are plotted on the horizontal axis and sequence 
data values are plotted in the vertical axis (see appendix). There exists a thresh- 
old value in r = 3.57 above which the transition of the dynamical system to the 
chaotic behavior can be observed. The bifurcation diagram serves as a reference 
to control results, and it is usually apphed in this resear ch frame 

Figure 5 shows results for the largest Lyapounov exponent estimate Ai, calcu- 
lated for each of the time series in the set, as a function of the control parameter 
r. This reference is taken as a control for the discussion of the results obtained 
by applying the DFT representation, and it is not a particular aim of this work. 

Figure 6 shows a mesh plot of the DFT representation amplitude spectrum 
abs{X[k\), (in order to obtain a finest version of each figure and figure showing 
a mesh plot of the DFT representation phase spectrum angle{X[k]), please, see 
appendix), for the discrete-time series matrix. The x axis corresponds to the 
control parameter r, and the y axis corresponds to the frequency, from zero 
value to the Nyquist frequency value. In order to obtain the best view of the 
DFT representation amplitude response mesh plot, the dc component on the 
signal was eliminated. 

3 Discussion of the results 

The application of the DFT representation of a chaotic-type signal gives satis- 
factory results for determining the order-chaos transition in the logistic equation 
system. This conclusion is confirmed by comparing that plot given in Figure 
6 with that shown in Figure 5, where the behavior of the largest Lyapounov 
exponent Ai with the control parameter r is shown. The largest Lyapounov 
exponent is one of the metric tools used for evaluating the system dynamics. 
Besides, bifurcation diagram corresponding to the logistic equation system may 
be used too as a reference. A critical value for control parameter r is reported 
in literature above which sequences produced by the logistic equation exhibit, 
all of them, a chaotic behavior. This threshold value is around r 3.57, as is 
shown in Figure 5. 

Some intermediate ordered region in the chaotic behavioral region underlie, 
but some authors [pT| refer them as a consequence of the computational system 
limitations rather than a real behavior of a dynamical system. However, this 
relationship between order and chaos can be observed from the dependence of 
the largest Lyapounov exponent Ai with control parameter, r, whenever the 
largest Lyapounov exponent value tends to zero at bifurcation points embedded 
in chaotic region (r > 3.57). The behavior of a dynamical system, for control 
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parameter values correponding to largest Lyapounov exponent close to zero, is 
not chaotic. 

The reconstruction parameter values selected are in correpondence with lit- 
erature, and satisfy requirements of the Takens theory and Grassberger and 
Procaccia algoritm[^. The bifurcation diagram was builded up several times, 
changing the initial condition, and changing the step value for the control pa- 
rameter, and the feature of the bifurcation diagram for r > 3.57 remains. It 
can be remarked that the order-chaos transition for the value r ^ 3.57 can be 
observed in DFT representation amplitude spectrum mesh plot (and in DFT 
representation phase spectrum mesh plot, see appendix), shown in Figure 6. 

On the other hand, Figures from 1 to 3 show the time domain representation 
of some of time series and their corresponding histogram. Comparing Figure 2 
(r = 3.70) and 3 {r — 3.83), it can be deduced the bifurcation effect from their 
histogram. Note that for figure no. 3 the control parameter satisfies r > 3.57, 
and it corresponds to a value for the largest Lyapounov exponent close to zero 
for r = 3.83, as it can be estimated in Figure 5. It can be remarked that the 
bifurcation diagram analysis is used in studies of the behavior of a dynamical 
system that tends to chaos through period-doubling procedure 

Applying DFT representation to discrete-time series in a family derived from 
a observable in a dynamical system is recommended for qualitative detecting of 
order-chaos transition, and this can be add to the metric tools of the nonlinear 
dynamical analysis. Any method that allow us to evaluate the behavior of a 
dynamical system, and that can be applied to a relatively short discrete-time 
series, must be taken into account as a practical method. 

4 Conclusions 

The application of the time series DFT representation can be made over a 
chaotic-type time series, nevertheless this time series is to be a nonlinear one. 
This method have not strong requirement about the length of the experimental 
data sequence. If one deals with a nonlinear dynamical system that can modify 
its behaviour between order and chaos, a time series produced by the system 
reflects this change, and the DFT representation of that time series exhibits in 
its amplitude response and its phase response the change in the system too. 
Any method, that allow us to give some evaluation about the behaviour of a 
dynamical system and that can work over a time series of a short length, adquires 
a practical importance. The DFT representation can be used for complementing 
the results obtained by applying the most powerfuU methods of the nonlinear 
dynamical analyisis. 
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5 Appendix 



In order to see the bifurcation diagram, a finest version of all figures including 
the DFT representation phase spectrum mesh plot, the following program may 
be runned on MatLab. For running on MatLab, write a percent symbol at the 
begining of each comment line. Comment: [Figures for DFTQI paper.] Com- 
ment: [To producing a family consisting of 120 time series, each of one having 
1024 data points.] close all, clear. Comment: [Datos, r = 2.8; del = 0.01; 
nr = 120; x = 0.65; iV = 1024;] for j = 1 : nr;R{j) = r;r = r + de?;end; 
for J = 1 : nr;fori = 1 : N;A{i,j) = x;x — R{j) * a; * (1 — a;);cnd;cnd; Com- 
ment: [Figures selected, j = [14 18 21 24 61 72 91 104 109];] for fc = 1 : 
9;rf{k) = i?(j(fc));end;rf, disp('values of j equals to 14, 94, y 104 were selected'); 
R = [2.93 3.70 3.83]; pause(l); figure(l); subplot(211); plot(A(l:128,14)); ti- 
tle('Figure l:Time scries and its histogram, for r = 2.93'); xlabcl('Timc index'); 
ylabelC Sequence values'); subplot(212); hist(A(l:128,14),40); xlabel('Sequence 
values'); ylabel('Relative frequency'); pause(l); figure(2); subplot(211); plot(A(l:128,94)); 
title('Figurc 2: Time series and its histogram, for r 3.70'); xlabel('Time in- 
dex'); ylabel('Sequence values'); subplot(212); hist(A(l:128,94),40); xlabel('Sequence 
values'); ylabel('Relative frequency'); pause(l); figure(3); subplot(211); plot(A(l:128,104)); 
title('Figure 3: Time series and its histogram, for r = 3.83'); xlabcl('Time in- 
dex'); ylabel('Sequence values'); subplot(212); hist(A(l:128,104),40); xlabel('Sequence 
values'); ylabel('Relative frequency'); Comment: [For optimal filling criterium.j 
r = 2.8; nx = 2^"; del = 0.01; nr = 120; for j = 1 : nr; R{j) = r; 
xO = 0.65; vO = 0.01; x = xO; v = vO; ioi i = 1 : nx; Al{i,j) = x; 
Bl{i,j) = v; z = x; x = r * x * (1 — x); v = x — z; end; r = r + del; end; 
pause(l); figure(4); subplot(221); plot(Al(l:nx,114),Bl(l:nx,114)); title('Figure 
4a : r = 3.94'); Comment: [Part of the program for largest Lyapounov exponent 
ploting.] xr = 2.82 : 0.02 : 2.8 + 60* .02; 



yexl = 


[-0.572 


- 0.240 


- 0.220 


-0.180 


- 0.153 


-0.118 


-0.091. 


-0.052 


- 0.031 


- 0.000 


- 0.055 


-0.127 


-0.201 


- 0.284 


- 0.381. 


-0.492 


- 0.625 


- 0.781 


- 0.996 


- 1.310 


- 1.250 


- 0.855 


- 1.030. 


-1.166 
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0.577 0.589 0.651 0.608 - 0.067 0.478 0.661 0.710... 

0.752 0.797 0.825 0.850 0.945] 



pause(l); figure(5); plot(xr,yexl,'kx'); hold on; plot(xr,yexl,'k:'); grid on; xla- 
bel('Control parameter'); ylabel('Largest Lyapounov exponent'); title('Figure 
5: Behavior of LLE with r'); legend( 'Experimental values', 'Behavior', O);zoom 
on; Comment:Part of the program for calculating and ploting the amplitude 
spectrum mesh and the phase spectrum mesh of the DFT representation. A = 
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A - ones{N, 1) * mean{A); f = linspace{0,pi,N/2); ftA = fft{A); amftA = 
absiftA{l : N/2,:)y, anftA = angle{ftA{l : iV/2,:)); pause(l); figure(6); 
mcsh(R,f,aniftA); titlc('Figurc 6: Amplitude spectrum mesh'); xlabel('Control 
parameter'); ylabel('Frequency'); zlabel('Amplitude'); view(-38,22); pause(l); 
figure(7); mesh(R,f, anftA); title('Figure 7: Phase spectrum mesh');xlabel('Control 
parameter'); ylabel('Prequency');zlabel('Phase');view(-40,30); 
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Figure 1: Time series and its histogram, for r=2.93 



Figure 2: Time series and its histogram, for r=3.70 



Figure 3; Time series and its histogram, for r=3.83 



Figure 4: Optimal filling of the phase space criterium 



Figure 5: Behavior of the largest Lyapounov exponent 



Figure 6: Amplitude spectrum mesh 
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Figure 1 : Time series and its Inistogram, for r = 2.93 
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Figure 2: Time series and its Inistogram, for r = 3.70 
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Figure 4 : r=3.94 
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Figure 5: Behavior of LLE with r 
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